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Computation of Transonic Flow over Porous Airfoils 

Thesis directed by Professor Chuen-Yen Chow 

A numerical tool is constructed to examine the effects of a porous 
surface on transonic airfoil performance and to help understand the flow struc- 
ture of passive shockwave/ boundary layer interactions. The porous region is 
located near the shock with a cavity underneath it. This study is composed 
of two parts. Solved in the first part, with an inviscid-flow approach, is the 
transonic full-potential equation associated with transpiration boundary con- 
ditions which are obtained from porosity modeling. The numerical results of 
this part indicate that a porous airfoil has a wave drag lower than that of a 
solid airfoil. The observed lambda-shock structure in the wind-tunnel testing 
can be predicted. Furthermore, the lift could be increased with an appropriate 
porosity distribution. 

In the second part of this work, the modified version of either an inter- 
active boundary layer (IBL) algorithm or a thin-layer Navier-Stokes (TLNS) 
algorithm is used to study the outer flow, while a stream-function formulation 
is used to model the inner flow in the shallow cavity. The coupling proce- 
dure at the porous surface is based on Darcy’s law and the assumption of a 
constant total pressure in the cavity. In addition, a modified Baldwin-Lomax 
turbulence model is used to describe the transpired turbulent boundary layer 
in the TLNS approach, while the Cebeci turbulence model is used in the IBL 
approach. According to the present analysis, a porous surface can reduce the 
wave drag appreciably, but can also increase the viscous losses. As has been 
observed experimentally, the numerical results indicate that the total drag is 
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reduced at higher Mach numbers and increased at lower Mach numbers when 
the angles of attack are small. Furthermore, the streamline pattern of passive 
shock/boundary layer interaction are revealed in this study. 
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CHAPTER I 


INTRODUCTION 


Background 


A method often adopted to reduce the drag associated with shock- 
waves and to improve performance envelopes for transonic aircraft is the use 
of supercritical shock free airfoil sections for the wings. Three computational 
•procedures have commonly been applied to design a supercritical airfoil: 

1) Procedures involving indirect methods. The hodograph and fictitious gas 
methods are in this category. 

2) Procedures involving inverse methods. Methods for solving the classical 
inverse problem of aerodynamics are in this category. 

3) Direct methods. This category is characterized by use of a direct compu- 
tational fluid dynamics (CFD) analysis program, coupled with a numerical 
optimization algorithm. 

For an extensive list of references and survey articles the reader is referred to 
Holst et al . 1 

However, supercritical airfoils are only effective at quite restricted 
design conditions, since drag increases rapidly at off-design conditions. 2 In 
order to extend the optimum conditions, a simple and economical concept for 
drag reduction was suggested by Dennis Bushnell and Richard Whitcomb of 
NASA Langley Research Center in 1979, according to which a passive shock- 
wave/boundary layer control is achieved by using a porous surface with a 
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plenum chamber underneath the shock location. This device can be catego- 
rized into combined suction and blowing devices, which have been used for 
boundary layer control (BLC) since 1940. 3 

Theoretically, the drag reduction of a transonic airfoil can be achieved 
if boundary layer separation and shock wave strength can be controlled by 
applying appropriate biowing or suction at the airfoil surface. Appropriate 
blowing in the supersonic region ahead of a strong shock may cause it to 
degenerate into a series of weaker waves or to generate another oblique shock 
upstream of the injection region, thus resulting in a smaller pressure gradient 
and a smaller entropy change. The additional kinetic energy supplied by 
blowing also increases the mixing rate in the boundary layer and prevents 
flow separation. However, strong blowing not only thickens the boundary- 
layer but also probably provokes an early separation as a side effect. On the 
other hand, the application of suction in the strong adverse-pressure gradient 
region would possibly delay separation but might produce a stronger shock 
and cause a higher wave drag as a side effect. In addition, if the suction area 
is of limited extent, it is necessary to examine whether the resulting boundary 
layer is capable of overcoming the adverse pressure gradient downstream of 
the suction region. Furthermore, according to an approximate nonasymptotic 
approach for weak shock/turbulent boundary layer interaction (including the 
mass transfer effects for attached flow), a small amount of suction could hasten 
the onset of separation slightly behind the shock foot. 4 Either blowing or 
suction requires power, thus an extra pump drag should be added to the total 
drag of the airfoil when an active control device is used. A passive control 
device (as sketched in Fig. 1.1), which provides blowing and suction without 
externally supplied power, hopefully, reaps the benefits of both blowing and 




Fig. 1.1 Porous airfoil in transonic flow. 
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suction without their negative effects and without excessive use of space. This 
expectation has been initially proven by the experimental results obtained by 
Bahi et al. 5 and Nagamatsu et al. 6-7 which indicated that at high-subsonic 
Mach numbers a properly arranged porosity can reduce total drag appreciably. 
But at lower Mach numbers, the drag was increased. Also a lambda-shock 
structure was observed above the porous airfoil in the experiments. In these 
experiments, the airfoil model was mounted so that the lower half of the airfoil 
was embedded in the bottom wall of the test section and only the upper surface 
was exposed to the transonic flow. The porous surface extended over about 
27% of the chord length. 

Savu et al. 8 presented Schlieren photographs for two airfoil models. 
The porous surface extended from near the leading edge to the trailing edge. 
They observed a reduction in shock wave strength as well. 

The passive shock wave/boundary layer interaction was investigated 
by Krogmann et al. 9 and Thiede et al. 10 They introduced a double-slot model 
and also studied the active suction effects. Their experimental results for a 
VFW VA-2 supercritical airfoil indicated that passive effects were effective 
primarily at relatively high Mach numbers and high angles of attack. In addi- 
tion, they reported both a reduction in the drag and an increase in the lift on 
a porous airfoil model. Furthermore, their results indicate that local suction 
in the shock region delayed the shock-induced separation and considerably 
improved the overall aerodynamic performance. For a double slotted config- 
uration at high angles of attack, no severe buffeting was noticed even though 
the flow was separated near the trailing edge. In general, they did not observe 
that the shock strength was significiently affected. These authors expressed 
uncertainty in the accuracy of their boundary layer probe measurements. 



Recently, Raghunathan and Mabey 11 studied a 6 percent thick cir- 
cular arc half airfoil set on a wind tunnel roof. They investigated the effect of 
hole orientation on the passive shockwave/boundary layer control and found 
that forward facing holes can reduce the drag appreciably. They also obtained 
a reduction in drag at higher Mach numbers and an increase in drag at lower 
Mach numbers for a porous airfoil, which agrees with the results of Bahi et 
al. 5 and Nagamatsu et al. 6 *' 

Before numerical studies of porous airfoils are surveyed, it is helpful 
to briefly review previous studies of transonic flow over solid airfoils. Although 
some physical assumptions may not be strictly valid and numerical instability 
may result from the treatment of boundary conditions, in principal, the nu- 
merical approach used for solid airfoils can be extended to the study of porous 
airfoils. 

The problem of transonic flow over a solid airfoil has been extensively 
studied with three different general approaches. The first approach is solving 
the nonlinear inviscid governing equations which can be either the transonic 
small-disturbance potential equation, the full-potential equation or the Euler 
equations. Only a small amount of computation time is required. This ap- 
proach provides a good preliminary computational tool that is valid in the 
absence of strong viscous-inviscid flow interaction. 

The second approach utilizes an interactive boundary layer (IBL) 
method. The essential components of the IBL method are an inviscid-flow 
algorithm coupled with a boundary-layer algorithm. The basis for this proce- 
dure was originated by Prantl, 12 who considered the high Reynolds number 
flow field to be mainly inviscid except for a viscous boundary layer near the 
airfoil and a viscous wake behind the airfoil. The principal interaction between 
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the boundary layer, wake and external inviscid flow arises from displacement 
thickness effects leading to a a thickened semi-infinite equivalent body which 
can cause significant changes in the surface pressure and forces. Separate 
systems of equations can be solved for the external inviscid flow and for the 
boundary layer and wake. The boundary-layer equations for the boundary 
layer and wake can be solved w’ith finite-difference or integral methods. The 
external inviscid flow can be represented by one of the nonlinear inviscid ap- 
proximations which have been mentioned earlier. The matching can be imple- 
mented at one of three different surfaces: 1) the surface of the airfoil and wake 
centerline, 2) the equivalent displacement surface, 3) the edge of the boundary- 
layer. Reviews of matching procedures for the two separate systems can be 
found in references by Melnik, 13 LeBalleur, 14 and Lock and Firmin. 1 ’ 

The third approach is to solve the mass-averaged Navier-Stokes equa- 
tions or thin-layer Navier-Stokes (TLNS) equations. The TLNS equations are 
obtained by neglecting all streamwise derivatives of the viscous terms, conduc- 
tive heat-flux terms, and any term involving mixed derivatives from the full 
Navier-Stokes equations. These neglected terms generally cannot be resolved 
with the available grid resolution. Recent work by Visbal and Shang 16 clearly 
indicates that the TLNS approximation produces essentially the same results 
as the full Navier-Stokes equations. This third approach (Navier-Stokes) is 
more expensive but imposes less theoretical restrictions than the other ap- 
proaches. 

Until now, transonic porous airfoils or transonic airfoils with mass 
transfer have been studied numerically only by the first two approaches. Savu 
and Trifu 17 and Savu, Trifu and Dumitrescu 8 have computed the flow about 
a porous NACA 0012 airfoil and a NACA 64A205.5 airfoil using the first ap- 
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proach (solving the transonic small-perturbation potential equation). Their 
results show that the standing shock can be eliminated completely by choos- 
ing the appropriate distribution of porosity for a given Mach number. The 
porosity was distributed from near the leading edge to the trailing edge. 

Mildly separated transonic flows over porous airfoils have been stud- 
ied by Oiling 18 via the second approach. An integral method was applied 
to the boundary layer equations and a full-potential solver was used for the 
inviscid flow. A slight improvement in airfoil performance due to porosity was 
found in this study. 

Shock wave/laminar boundary layer interaction with mass transfer 
(via the second approach) has also been studied by Ram et al . 19 The boundary 
layer equations were solved by an integral method with attention to necessary 
modification in boundary conditions for flows with mass transfer through the 
surface. The zonal approach was adopted for solving different equations in 
different zones. The results indicate that full-chord laminar flow can be main- 
tained and separation can be prevented by the use of suitable suction. 

Recently, Inger and Nandanan 20 presented a nonasymptotic triple 
deck approach for shock/turbulent boundary layer interaction (SBLI) that 
included the effect of wall suction confined to the SBLI zone. It was found 
that both displacement thickness and momentum thickness were increased 
along the. aft portion of airfoil at low suction rates. This work (employing the 
second approach) was limited to unseparated flow, and suction effects on the 
shock position were neglected. 
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Motivation 

The aforementioned experimental and theoretical results indicate 
that porous surfaces can affect an airfoil performance at transonic speeds. Un- 
like a supercritical airfoil whose fixed shape is designed for certain optimum 
Mach number and angle-of-attack ranges, a porous airfoil can be adjusted 
without changing its contour shape. The adjustments achieved from variable 
surface porosity may result in improved performance over a wide range of 
conditions. Experimental demonstration of the possible gains has been ham- 
pered by uncertainties arising from wall and support interference in transonic 
Wind tunnels. Theoretical evaluation is impaired by the absence of a complete 
theory for flow over a transonic airfoil with mass transfer. In view of the 
importance of the problem and the improving computer capabilities, it seems 
worthwhile to continue the numerical investigations. 

The purpose of the present study is to achieve improved computa- 
tions of transonic flow over porous airfoils. An inviscid code, an IBL procedure, 
and a TLNS procedure are employed. To the author’s knowledge, no previous 
work has brought these procedures to bear on the problem at hand. 

The inviscid code TAIR, originally developed by Holst, 21 solves the 
conservative, transonic full-potential equation using an approximate factoriza- 
tion algorithm (AF2). The IBL procedure TRIVIA, developed by Van Dalsem 
and Steger, 22-23 consists of a direct-inverse, finite-difference boundary-layer 
algorithm coupled with TAIR. The TLNS procedure ARC2D, developed at 
NASA Ames Research Center, 24 is based on the Beam and Warming implicit 
approximate factorization algorithm with several further improvements. 

Chapter II includes a description of the inviscid governing equations, 
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solution algorithm, and treatment of boundary conditions for transpiration 
at the airfoil surface. In Chapter III, the viscous boundary layer equations, 
the modified algorithm for the boundary layer equations, the modified viscous- 
inviscid interaction algorithm, and the Cebeci turbulence model are described. 
Chapter IV includes a description of the TLNS equations, the algorithm, and 
a modified version of the Baldwin-Lomax turbulence model. The cavity-flow 
solver and porosity modeling are described in Chapter V. The numerical re- 
sults are discussed in Chapter VI. The conclusions from this work and recom- 
mendations for future work are presented in Chapter VII. 



CHAPTER II 


INV1SCID FLOW SOLUTION 
Transonic Full-Potential Equation 


Recause the changes of density, velocity and temperature across a 
shock are first order but the entropy jump across the shock is third order, 25 
inviscid transonic flow can be reasonably approximated as an isentropic and 
irrotational flow when the Mach number ahead of the shock (M\) is less than 
1.3. 26 This assumption allows the flow to be described by a potential <J> which 
satisfies the continuity equation 

{p$x)x + [p$y)y = 0 (2.1a) 


7 ~ 
7 + 




(2.16) 


where x and y are Cartesian coordinates nondimensionalized by the airfoil 
chord length c; p and 3> z ,$ y are the density and velocity components nondi- 
mensionalized by the stagnation density p? and the critical sound speed a' , 
respectively; and qf is the ratio of specific heats. 


For computational convenience, the governing equations are trans- 
formed from the physical domain (Cartesian coordinates) to a computational 
domain by the transformation 



1 ] 


£ = £(z,y) v = v{x,y) 

The full-potential equation written in the £-77 computational domain is given 

by 


where 



7- 1 
7 + 




(2.2a) 


( 2 . 26 ) 


U — A]$c + A2^rj V" — _i4 2 $£ + 


^t = e* + ej 

■^2 T ^y^ly 

^3 = 7x + 

and 

The variables £/ and K are the contravariant, velocity components along the 
£ and 77 directions, respectively; A j , ,4 2 , and ^3 are metric quantities; and 
J is the Jacobian of the transformation. This transformation maintains the 
strong conservation-law form of the original equation and hence possesses 
characteristics suitable for a shock-capture scheme. The Ai and As metrics 
provide a measure of cell aspect ratio. The A\ metric is approximately the 



12 


ratio of arc length along the ^-direction to arc length along the ^-direction. 
The A 2, metric is approximately the inverse of A\, and A2 is a measure of 
orthogonality. The .Jocobian. .1. can be shown to approximate the inverse of 
cell area. 


Numerical Algorithm 


The numerical algorithm and code. TAIR, developed by Holst 21 are 
capable of simulating inviscid flow about an airfoil (with weak shocks) by 
solving the transonic, full-potential equation in body-fitted coordinates using 
an AF2 scheme. This code provides rapid convergence and requires only a few 
seconds of computer time per case on a CRAY-XMP processor. 

The AF2 fully implicit scheme can be expressed as 

(a - My)(a*; - JcAjc)^ = auL^ (2.3) 


in which C^ - = $” J hl - ft is an acceleration parameter, w is a relaxation 
parameter (2 > w > 0) and A, and Aj are defined by 

A < = (t 1 ) . 

The scheme is implemented in a two-step format. In step 1, a scalar bidiagonal 
matrix is inverted for each £ — constant line and in step 2, a scalar tridiagonal 
matrix is solved for each rj = constant line. 

The residual on the right-hand side of Eq. (2.3) is given by 

7 d P J). , . + *« ,(-,-) , 

J l +5J J 


where 6 and 6 n are first-order accurate, backward-difference operators in 
the f and r\ directions, respectively. The artificial density scheme consists of 
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introducing an upwind evaluation of the density in f and tj directions (denoted 
by p and p). Basically, the supersonic regions are stabilized using an upwind 
bias of the density. This provides an efficient and reliable spatial differencing 
scheme for the capture of weak (transonic) shock waves. Details of this scheme 
can be found in Refs. [21,27]. 

There is an entropy-correction option in the TAIR code which in- 
volves replacing the isentropic density p t by 


P = 


Pie 


-AS/H 


The normalized entropy change (A S/R) can be approximated by a locally 
one-dimensional shock relation which is a function of M\ (the Mach number 
upstream of the shock). 


AS 

~w 


2 1 


3 ( 7 + 1 ) 


2 



(2.4) 


Another'feature of the TAIR code is its ability to capture free-stream 
flow in general curvilinear coordinates. Three free-stream consistency condi- 
tions are required for free-stream capturing. The reader is referred to Ref. 
[28 j for complete details. 

The body-fitted grids used for the full-potential solver are generated 
by the finite-difference solution of Poisson equations using a computer code 
(GRAPE) developed by Sorenson 2 '* and based on the work of Steger and 
Sorenson. 30 Two sets of Poisson equations are solved using the successive line 
over-relaxation procedure. This grid generation code allows control of the 
grid point spacing along and normal to the boundaries as well as the angles 
at which grid lines intersect the boundaries. The C-mesh grid topology used 
in this study is shown in Fig. 2.1. A typical inviscid grid is shown in Fig. 2.2. 
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Fig. 2.2 Grid used for inviscid-fiow calculation. 
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Boundary Conditions 

The boundaries associated with the physical domain are transformed 
to boundaries of the computational domain. At the outer boundary, the ve- 
locity potential is set to be the sum of the uniform free-stream component and 
the component caused bv a vortex with circulation T, where T is the jump 
in the velocity potential at the airfoil’s trailing edge. At the airfoil surface, 
there are two types of boundary conditions depending on the surface property 
(solid or porous). 

1. On solid regions of the airfoil, the flow tangencv condition is 
satisfied by requiring that the contravariant velocity component in the re- 
direction vanish. 

V = 0 on solid surfaces (2.5) 

2. On porous regions of the airfoil, the normal wall velocity v„ does not. vanish 
and its value is determined by porosity modeling, which will be introduced in 
Chapter V. 

The physical transpiration velocity v n is then transformed into the 
computational domain: 


v~ 


\fAz 

0 

Vn 

u 


Aa 

L ^ 

J 

J 

0 


where the negative sign is due to the use of a left handed coordinate system. 
With the C-mesh topology, the 77 -coordinate lines intersect the body at close 
to right angles and, to a good approximation, A 2 — 0 at the body. Using 
this approximation v n does not contribute to U and can be expressed in the 
computational domain as: 


V = -\f~AzVn 


on the porous surface 


(2.7) 
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Numerical implementation of the boundary conditions on the airfoil 
surface is summarized as follows: 

1. On the solid portion of the airfoil, the flow tangency boundary 
condition (i.e., V = 0) is used to obtain 


J ) i,NJ + i \ J / i,N J - 

where j ~ N J is the airfoil surface (Fig. 2.1). 


( 2 . 8 ) 


2. On the porous portion of the airfoil, the boundary condition (i.e., 
v - -VT 3 v n ) is applied with the aid of a Taylor series expansion 


PL 

J 


pV 


± 


A ij ( pV 

T 


i.NJ 


i,NJ ± i \ ^ / i,NJ 

+ ± 0|(Ai)) 3 ] 
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J 


i,NJ 


so that the first-order boundary condition becomes 


(2.9) 


pV 


pV 


J ) i,NJ + k V J J i,NJ --i V J / i,NJ 


pV 


( 2 . 10 ) 


In fact, Eq. (2.8) is a special case of Eq. (2.10) for V — 0 at j = N J . 

Equation (2.2), together with the surface boundary conditions [Eqs. 
(2.8) and (2.10)1 and appropriate outer boundary conditions, are solved nu- 
merically using the modified TAIR computer code. At the beginning of each 
iteration, the porous boundary condition is updated by porosity modeling. 
Once the flow field is obtained, forces on the airfoil in the x and y directions 
can be computed by performing integrations around the airfoil contour 


= / [pn x + pv n u)ds 

Js 


( 2 . 18 ) 
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in which S is the circumference of the airfoil; u and v are velocity components 
in the x and y directions, respectively; p is pressure; and n x and n v are direc- 
tion cosines between the airfoil surface and the x and y directions, respectively. 
Equations (2.18) and (2.19) are derived from the momentum conservation the- 
orem. Note that the second term in both equations ( pv n u and pv n v) are zero 
for solid wall airfoils but must be retained in the present calculations because 
of the porous wall assumption. Once F x and F y are obtained from Eqs. (2.18) 
and (2.19), the lift and drag coefficients, Ci and Cp, are easily computed. 

The drag components of a two-dimensional airfoil are wave drag, 
viscous pressure drag, friction drag and pump drag. At this stage, the only 
drag component evaluated is wave drag because viscous effects are neglected 
and no power is supplied. 


/.< 


( pn y + pv n v)ds 


(2.19) 


L 



CHAPTER 111 


INTERACTIVE BOUNDARY LAYER PROCEDURE 


A finite-difference viscous-inviscid interaction algorithm is modified 
to simulate the viscous transonic flow over porous airfoils in the presence of 
mild separation. The code consists of a direct-inverse, finite-difference bound- 
ary layer algorithm coupled with the full-potential solver described in Chapter 
II. A viscous-inviscid interaction algorithm is used. For the present work, ma- 
jor modifications in this procedure are associated with the forcing function, 
the boundary condition treatment, the viscous/inviscid interaction algorithm, 
integration of the continuity equation, and the turbulence model. 


Governing Equations 

The nondimensionalized, first-order compressible turbulent bound- 
ary layer equations for the steady, two-dimensional flow of a perfect gas are: 


x-momentum equation 

p(uu x + UUy) = -f3p x + ((pi + fJL t )u y )y (3.1a) 

energy equation 

pc p (uT x + vT y ) = fiup x + ((«/ + k t)T y )y + {m + p,t){uy) 2 (3.16) 

perfect gas equation 

P — pT (3.1c) 

continuity equation 

(pu) x + [pv) y 0 (3. let) 
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where x,y are along and normal to the body or wake centerline, respectively. 
Viscosity, pressure, temperature, and density are nondimensionalized by their 
free-stream values. The u and x variables are nondimensionalized by the free- 
strcam velocity and a characteristic length, respectively, whereas the v and y 
variables are nondimensionalized by the same quantities divided by \/Re oq. In 
Eq. (3.1), (3 = (p/pU 2 ) oo, pi and k. / are molecular viscosity and conductivity 
determined by the Sutherland viscosity law and a constant Prandtl number 
assumption, pt and /c* are eddy viscosity and conductivity determined bv tur- 
bulence modeling. These equations are simplified from the Reynolds equations 
with use of the Boussinesq approximation, that is, 

-pu'v' — p t u y (3. 1 e) 

-Cppv'T' = K. t T y (3.1 /) 

Experiments confirm that the ratio of the diffusivities for the turbulent trans- 
port of heat and' momentum, called the turbulent Prandtl number, Pr t - 
p t c v /K.f is a well-behaved function across the flow. Most algebraic turbulence 
models set Prt — 0.9 such that only pt needs modeling. 

To solve the boundary-layer equations using finite-difference approx- 
imations, it is necessary to construct a grid. In the interest of accuracy and 
computational efficiency, it is desirable to use a non-uniform grid. However, if 
Eqs. (3. la-3. Id) are differenced on a non-uniform grid, complicated variable- 
spacing finite-difference operators must be used. On the other hand, if a gen- 
eral x, y to £(x),y(x,y) coordinate transformation is applied to the boundary- 
layer equations, a non-uniform x,y grid may be used while the equations can 
be solved on a uniform grid. In other words, the physical domain x,y grid 
is adapted to resolve the flow field, while the computational domain grid 
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is chosen so that simple equal-spaced finite-difference operators may be used. 
Hence, although the equations are slightly more complicated after the coor- 
dinate transformation, overall the finite-difference solution of the equations 
becomes much simpler. 

For computational convenience, the boundary-layer equations are 
transformed from the physical x,y domain to a uniform £,(x),r](x,y) com- 
putational domain: 

x-momentum equation 

P\u{uc£ x + u v ri x ) + vu v r] y \ = -/3pc£ x + ((in + n t )u rt n y ) T ,r)y (3.2a) 
energy equation 

Pc p \u( T s( x + Tr,Ti x )+vT n riy\ = 0upc£ x + ((K l + K t )T T ,r]y) r) 7i y + (ni + Ht)(u„Vy) 2 

(3.26) 

perfect gas equation 

p=pT (3.2c) 

continuity equation 


{pu)i£ x + (pu) n T] x + (pv) v 'n y - 0 (3. 2d) 

Moreover, it is also frequently assumed that the stagnation enthalpy 
is constant across the boundary layer in the case of low- to moderate- speed air 
flow over an adiabatic surface. This approximation does not introduce large 
errors provided the Mach number M e at the boundary-layer edge is moderate 
(the actual variation of stagnation enthalpy across an adiabatic boundary- 
layer is approximately equal to 4% when M e = 2). Therefore the energy 
equation can be replaced by a simple algebraic relation ( H = c p T + u 2 / 2 = 
constant), and an increase in computational efficiency achieved. 
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Numerical Algorithm 

The finite-difference boundary-layer equations are solved using a 
predictor-corrector space marching algorithm with appropriate initial condi- 
tions as shown in Refs. j22-23). However, integration of the continuity equation 
from the wall to obtain v ( v n can be nonzero) has been replaced by 

) + rix6 n {P u ) 

Vy 

where E is the shift operator (e.g. E n 1 u/ c = Uk-,), <5^ =1 — E'~ l , 6c = 
{Ec - E7 1 )/ 2, and 6 V = ( E ^ — E~ l )/2 . This procedure is simpler than the 
.original version. Since the boundary-layer equations are weakly coupled, they 
are solved in a sequential manner at each streamwise station. Overall, second- 
order accurate solutions are obtained at the cost of a few scalar bidiagonal 
and scalar tridiagonal matrix inversions per streamwise station rather than 
at least one block tridiagonal inversion per marching step of the box scheme. 
Both direct and semi-inverse interactions are built into the code. For attached 
flow, pressure is specified in the direct mode. However, near and in the reverse- 
flow regions, in order to avoid the Goldstein singularity, the wall shear t w and 
wake-centerline velocity u wc are specified in the inverse mode. 

By applying the x-momentum equation at the porous wall (u w = 
0 ,Un 7 ^ 0 ), one can obtain a relation : 

P^xP$ {pV-qTjy ) r^TJ-y |m pWUylly |u, (3-4) 

This expression allows us to eliminate the term /3(, X P$ from the differ- 
enced x-momentum equation and to put t w into that equation by the following 
approximation: 


(3.3) 


fir, ( P V ) = - 


1 + E n 1 Zxfi‘{pu 
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where subscripts 1 and 2 correspond to the wall and one grid away from the 
wall, respectively. Since these inverse forcing functions are cast in a general 
form, they can be applied to either a solid surface or a porous surface. In 
the inverse mode, t w is updated by the following viscous-inviscid interaction 
algorithm : 


C 1 (3-6a) 

-The acceleration parameter u is gradually increased to avoid the high fre- 
quency error which occurs during the first few iterations. A similar procedure 
is used to update the wake centerline velocity. 


c^c+^-prx 


(3.66) 


The quantity varies in a range near ten and oj c £ x varies near two. 

Viscous effects are introduced into the inviscid solver via a transpi- 
ration velocity v n determined from the boundary-layer solution 


1 / d(p e u e 8) 


Pu>Vr 


(3.7) 


pc \ ds 

where v n is the surface blowing (or suction) velocity determined from porosity 
effects, p e and u e are the inviscid values at the airfoil (or wake centerline). 
This equation is derived by considering the difference between the continuity 
equations for inviscid and viscous flow, 15 that is 


d 3 . Sr, 

— (PtU, - pu) + — ( piVi - pv) = 0 


(3.8) 
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where subscript i denotes inviscid quantity. Integrating across the boundary 
layer, from n = 0 to 6, one can obtain Eq. (3.7). since both (p,t>, — pv) and 
[piUi — pu) vanish at n = 6 due to the matching condition. The transpira- 
tion velocity v n is imposed at the airfoil surface in the inviscid flow solver 
(see Fig. 3.1). This interaction procedure avoids supercritical behavior and 
the the need for inviscid-grid generation at each iteration. As LeBallcur 31 
indicated, the supercritical behavior has no physical significance and can be 
controlled solely by the choice of the matching condition coupling the inviscid 
and boundary layer equations. In Fig. 3.1, the f l , r] 1 coordinate system is for 
the inviscid-flow solver and the £, r\ solution-adaptive coordinate system is for 
the boundary-layer solver. The total grid height for the viscous turbulent flow 
is a function of computed displacement thickness: 

ymax = oS'i-l ( 3 - 9 ) 

where O rv O . The viscous grid is generated by an exponential stretching 
function: 

Vj — Vj-i + Ay m m(l 4 t) 1 2 (3.10) 

The first grid point above the airfoil is placed at approximately y~ -~ 1 (which 
determines t). 

Once v n is known, the numerical boundary condition is applied as 
in Chapter II except v n in Eqs. (2. 7-2. 8) is replaced by v n . The complete 
interaction scheme is indicated in Fig. 3.2. From this IBL procedure, the skin 
friction and viscous pressure drags are evaluated in addition to the the wave 
drag. 
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Turbulence Model 
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The turbulence model used in the IBL procedure is the Cebeci two- 
layer algebraic turbulence model. The Prandtl formulation used in the inner 
layer is, 

o du 

(*<<)„,„„ = (3.11) 

The mixing length , / , is evaluated from 

/ = /cy(l — e y ! A ) (3.12a) 

The quantity y~ is given by y + = ya r /i/, and the friction velocity u T by 
•u T — {t w /p) X/ ~. In Eq. (3.12a) K. is the von Harman constant usually taken 
as 0.41 and is the damping constant most commonly evaluated as 2G. 
Cebcci- extended Van Driest’s modeling of the viscous sublayer and let 
be a function of , p ^ , p e and p e (subscript e denotes the edge of boundary 
layer). Hereu + = v n /u T and p + = (dp / dx) (l/ / pu T 3 ) . The Clauser formulation 
used in the outer layer is 

[vdouter = 0.0168pu e <5' (3.126) 

where u e is the velocity at the edge of boundary layer and 8" is given by 

rf 

S' = / (1 - — )dy (3.12c) 

Jo U t 

The outer value (Ht) ou t eT is used at all values of y beyond the point where 
-j[ Ut)inner an d (^Oouter ^ rst cross - The reader who is interested in this model 
is referred to Ref. (32) for further details. 


CHAPTER IV 


THIN-LAYER N AVIER-STOKES PROCEDURE 

The thin-layer Navier-Stokes equations are generally referred to in 
the literature as TLNS equations. A main feature of these equations is the 
presence of a nonzero, normal-pressure gradient, which is necessary to cou- 
ple and solve simultaneously the viscous and inviscid regions. By comparison 
with the Navier-Stokes equations, these composite equations require less com- 
putational effort because they contain fewer terms. However, the thin layer 
approximation is invalid at low Reynolds numbers and in regions of massive 
flow separation. 


Governing Equations 

Again for computational convenience, the governing equations are 
transformed from Cartesian coordinates to general curvilinear coordinates by 
the transformation 

t — t 

f = f(x,y,t) (4.1) 

T] = T](x,y,t) 

The conservative thin layer Navier-Stokes equations expressed in terms of 
Reynolds mass-averaged variables in general curvilinear coordinates (£,77) and 
time r can be written as 


d T Q + dcE -|- d n F = Re 1 d n S 


(4.2) 
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where 

' P "I r pu 

n - r-i P u r - 7-1 puf/ + fxP 
pv pvU + ZyP 

. e UJ(e + p) - &p 

P V 

puV + r) x p 
PV V -f PyP 

V(e-r p) - T) t p 

with 

U = & + + fyV 

V ~ 711 + Vs.v. + rj y v 

and 

0 

+ r) y m 2 
r) x m 2 + r/ v m ; » 

7 } x (u m i + tim 2 + m 4 ) + p y {mn 2 + V 7 n :i + 7715) 

wherein * 

mi = p(47?x^rj - 2pyVr 1 )l'i 
m 2 = p{ri y Un + Px^) 

= p[~2p x Uy + 4p y v n )/3 
m 4 = pPr~ l (-j - l) _1 7/ z 5 7? (a 2 ) 

m 5 = /xFr -1 (7 - l) _ l p y dy(a 2 ) 

Pressure is related to the conservative flow variables, Q, by the equation of 
state 

P=(7 — 1) e - ^P(« 2 + ^ 2 ) ( 4 - 3 ) 

where 7 (=1.4) is the ratio of specific heats. The speed of sound is a (for ideal 
fluids, a 2 = 7 p/p). The dynamic viscosity p is made up of an eddy viscosity 
(p t ) and molecular viscosity ( pi ) which is evaluated by Sutherland’s semi- 
empirical formula. Re is Reynolds number and pPr~ l represents p\Pr~ l + 
PtPr~\ 
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Numerical Algorithm 


The TLNS procedure AR.C2D developed at NASA Ames Research 
Center 15 is modified for this study. This implicit code was based on the Beam 
and Warming approximate-factored algorithm. Euler implicit or three-point, 
implicit time difference and second-order central spatial difference are utilized. 
An artificial dissipation model was added to capture a non-oscillating shock 
and maintain stability. Pulliam and Chaussee 33 developed a diagonal version 
to reduce the block pentadiagonal inversion to 4 x 4 matrix multiplications 
and scalar pentadiagonal inversions. Also a local time step and variable grid 
spacing are employed to accelerate the convergence rate for the steady How 
solutions. The reader who is interested in the numerical algorithm is referred 
to Rcf.'flS! for details. 

Boundary Conditions and Grid 

On the airfoil surface, the normal wall velocity v n is specified and 
the tangential velocity is set to zero. The boundary condition on pressure at 
the airfoil is taken to be dp/dr] = 0, since the rj coordinate lines are nearly 
orthogonal to the airfoil surface. The adiabatic wall condition is used to obtain 
density at the surface and total energy is decoded from the equation of state. 
The far field boundary condition is set by imposing a compressible potential 
vortex solution on the free stream quantities. 34 All the boundary conditions 
are updated explicitly in the TLNS procedure. 

The grid for the TLNS procedure is obtained by a hyperbolic grid- 
generation procedure. 35 Since the equation system used to generate the grid 
is hyperbolic in the tj direction, the outer boundary is not specified in advance 
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as in the elliptic grid-generation procedure. A typical grid set for the TLNS 
is shown in Fig. 4.1. By comparison with the inviscid grid (Fig. 3.2). the grid 
spacing for the TLNS procedure is more clustered near the airfoil surface. 
For most of the test cases in this work, the first point above the airfoil is at 
y + ~ 0(1), the total number of grid points is 251 x 65, and the outer boundary 
is set at 16 chord lengths away from the airfoil. 

Turbulence Model for Transpired Boundary 


The Baldwin-Lomax model' 30 used in the TLNS procedure is a two- 
layer algebraic model in which p t is given by 


[Pt) inner > 2 / — Ucroj.<ovcr 

(Vt) 0uler , y > y crossover 


(4.4) 


where v is the normal distance from the wall and y C rossover is the smallest 
value of y at which values from the inner and outer formulas are equal. 

The Prandtl-Van Driest formulation is used in the inner region 


Miner = (4.5) 


where 

l = ky\ 1 - exp ( -y + /A + )] (4.6) 

|w| is the magnitude of the vorticity and 

-f- PwU T y y/ Pw T wV 

V = — — = * (4.7) 

Pw Pw 

The eddy viscosity coefficient in the outer layer is given by 

{p l) outer = KCcppFw AKE^KLEB (y) (4-8) 
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Fig. 4.1 Grid used for TLNS calculation. 
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where K is the Clauser constant(0.0168) and Cep is an additional constant 
(1.6), and 


Fwake = min 


Vmax Fmaz 

('W KlhnuxUmF 2 /F„ 


(T9) 


^ DIF — M-i it y„,,, r Umin (4.10) 

where Cwk = 1-0, and y mtlT and F mrix are determined from the function 

F (y) - y\uj\ ll - exp(-y : /A + ) } (-1-11) 

The quantity F mai is the maximum value of F (y) that occurs in the profile, 
and i/max is the value of y at which it occurs. The Baldwin-Lomax model is 
patterned after that of Cebeci with modifications that avoid the necessity for 
(inding.the outer edge of the boundary layer. The length scales arc determined 
bv the distribution of vorticity. Since the original Baldwin-Lomax turbulence 
model did not consider blowing (or suction) efTects on A + , the damping factor 
A + is redefined as : 37 
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F+T 


(4.12a) 


where r + = r/r,„. In order to compare with the STAN-5 results in Refs. [38-39] 
near the blowing (or suction) surface, the x-momentum equation is integrated 
with v = v n = constant , p x = constant , p = constant and uu x neglected. 
The result is 


7-+ « t w + + p + y+ + v + u + {y + ) (4.12 b) 

where p + , y + , v + and u + are already defined in Chapter III. Fig. 4.2 shows 
plots of A + vs. p + with as parameter. The symbols represent data from 
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STAN-5. It is observed from Fig. 4.2 that when p~ « —0.04, and v + = 0, A~ 
goes to oc. For agreement with this, the denominator of Eq. (4.12a) should 
go to zero under those conditions 

r + = \ - 0.04 y + = 0 

Thereby y + = 25 is determined and extended to other cases with v + ^ 0. The 
results from setting y + = 25 in Eq. (4.12b) and n = 0.7 in Eq. (4.12a), [pt is 
evaluated with Eqs. (4. 5-4. 6)) are shown as the solid lines in Fig. 4.2. It is seen 
that values of A~ vs. for various u + from this procedure agree well with the 
results from STAN-5. With this modification, the dependence on boundary- 
layer edge quantities in the Cebeci turbulence model can be neglected. Indeed, 
this value of n is in the range proposed by other researchers, such as Patankar 
and Spalding (n = 0.5) 40 and Baker, Jonsson and Launder(n = 1.0) 41 . This 
modification (Eq. 4.12a) has been added to the Baldwin-Lomax turbulence 
model, yielding skin friction values in the blowing region higher than those 
calculated by the original Baldwin-Lomax model. Conversely, this model pre- 
dicts skin friction in the suction region lower than from the original model. 

For example, for the RAE2822 airfoil at M 0 0 = 0.7, a = 2.0 and 
Re = 6.5 x 10 G , with the transition position specified at 0.03 chord, and 
with the specified distribution of blowing shown in Fig. 4.3, the skin friction 
is reduced in the blowing region. But the modified version predicts slightly 
higher skin friction values than the Baldwin-Lomax model . The blowing 
effects on the C p distribution are plotted in Fig. 4.4. The effect of blowing 
on the pressure predicted by the two models is the same. Next the effect of a 
specified suction distribution is investigated and again the two models predict 
the same effect on the pressure, but the increased skin friction over the suction 
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region is slightly lower in the modified version shown (see Fig. 4.5). That is 
because in the modified version the blowing increases the mixing rate (.4' r 
being reduced) and the suction tends to iaminarize the (low (,4 + being in- 
creased). These effects are neglected when /l 4 ' = 26 is assumed. Similarly, a 
strong adverse pressure gradient would reduce A ~ according to the modified 
version. The skin friction C f\ e for the RAE2822. at Mcc = 0.73. Cl = 0.803. 
and Re = 6.5 x 10 G with transition specified at 0.03 chord, is shown in Fig. 4.6. 
Results from the two models are compared with Mehta’s' 12 computation and 
with experimental results 43 by Cook et al. Ail of the predictions agree well 
with experiment, but there is a slight improvement in results from the modified 
version. 

Surface roughness also has a large effect upon /l 4 . But this effect is 
not considered in the present work. 

Generally, the present IBL procedure is one to two orders of magni- 
tude faster than the TLNS procedure. For a mildly separated transonic flow, 
it can provide very good results. 22-23 However, the range of allowable blowing 
rates in the IBL procedure is more restricted than in the TLNS approach. This 
restriction arises mainly from the boundary layer theory assumption that the 
magnitude of the normal velocity component v is less than 0{Re~ x ^ 2 ). In ad- 
dition, the first-order boundary-layer approximation becomes suspect when a 
strong shock boundary-layer interaction occurs, because dp/dy in the bound- 
ary layer may not be negligible. At the expense of longer computing time, 
the TLNS procedure has less theoretical restrictions than the IBL procedure. 
Therefore, when a strong interaction occurs or blowing (or suction) becomes 
strong, the TLNS procedure is preferred in the present study. 
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Fig. 4.6 Comparison of C j\e distributions for an RAE 2822 airfoil 
at Mr* = 0.73, Re = 6.5 x 10 6 and C L = 0.803. 





CHAPTER V 


POROSITY AND CAVITY-FIOVV MODELING 

Since the holes over the porous region are very small and dense, it is 
difficult to compute the flow across the porous media without extremely fine 
grids. Therefore the porosity effect is modeled instead of being computed in 
the present study. The model proposed is patterned after the treatment of 
porous wind-tunnel walls ' 1 ' 1 based on the Darcy’s law. 

Interface of Outer and Inner Flows 


There are two porosity models used in this study. For both of them, 
the transpiration velocity v n for the outer flow is governed by Darcy’s law 
such that 


Vinner ) (5.l) 

where the subscript n indicates the direction of the outward normal on the 
surface; Vouter and p !rmer are the pressure above and below the porous surface, 
respectively; a is the porosity distribution function which is determined by 
viscosity as well as by the size and density of the holes in the porous surface. 
The subscript outer indicates the outer-flow property and inner indicates the 
inner-flow property. The first cavity model considered assumes a constant 
pressure, p, in the cavity, so that, 
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Vn — {.Pouter Pinner ) (5-2) 

Pc o'-' oc 

The constant cavity-pressure assumption was made for convenience but is 
expected to be qualitatively accurate. For passive flow through the cavity, the 
net mass flow through the porous surface of length 5 must be zero, or 


J (P^n) outerds — 0 


which leads to the relation. 


_ _ /$ °{PP) outerds ^ 

Pinner — r I ] 

J c OPouterdS 

Once pmner is known, the v„ for the outer flow can be determined from 
Eq. (5.2). 

For the second model considered, the total pressure in the cavity is 
assumed to be constant, that is, 


constant 


As before, total mass flow through the porous surface is conserved, and Darcy’s 
law is applied across the interface, 


fn — jj Pouter ( ) inner{Pt) inr 

Pco^ oo Pt 


which leads to the relation 


(P^) inner 


f s o{ P p) outerds 
0P O uter{ 'fr) innerds 


The pressure variation in the cavity [p/pt) inner can be determined by com- 
puting the cavity flow. For convenience in coupling with the outer flow, the 
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flow in the cavity is assumed isentropic and irrotational with a different total 
pressure from the outer flow. The total temperature is assumed to be the 
same as that of the free stream. Therefore the cavity flow can be described 
by the stream-function formulation 



(5.8a) 


P = 


7 — f (*s + *;) V 1, 

7 t1 p 2 


where 


pu = pv = -vp x 


(5.86) 


(5.8c) 


The stream-function is set to zero at both the vertical side walls and the 
bottom wall of the cavity. Along the porous surface, it is determined by the 
integration of mass flow rate from the outer flow by 

'f'(x) = f [pv n ) outer ds (5.9) 

J Xi 

as indicated in Fig. 5.1. 

The (p j p t ) inner > ( Pt)inner and v n can be determined iteratively. The 
iterative process uses the first model to obtain an initial v n , then integrates 
the mass flow rate by Eq. (5-9) to obtain the boundary conditions for the 
cavity and solves the cavity flow to update [p/ Pt)inner underneath the porous 
surface. Then update the {Pt) inner by Eq. (5.7) and updates v n by using 
Eq. (5.6) to finish the first iteration. Then iterative procedure is repeated to 
update the boundary stream function until the converged solutions (( Pt)inner , 
{p/pt) inner, an d v n ) are obtained. Then v n is available for the next update of 
the outer flow solver. 


i 



OUTER FLOW 
SOLVER 


DARCY'S LAW 



Fig. 5.1 Porous airfoil with a cavity-flow model. 
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The nonuniqueness of Neumann problem is avoided and boundary 
conditions are easily implemented by using the stream-function formulation 
(rather than full-potential equation with specified normal derivatives at all 
the boundary). The resulting model for the cavity flow, although not perfect, 
helps by providing estimates of how important pressure variations within the 
cavity might be. 

Cavitv-Flow Solver 

The stream-function formulation in the cavity can be written in the 
£-r) domain as 


where 


■hv - *iu 

r- j-' 


.4 3 .4 2 

T- b ~ T- x 


= o 


(5.10a) 


p = 


, 1 ~ 1 •; 

1 q 

7 + 1 


(5.106) 


U 




V = — — vi/ 
P 


(5.10c) 


./4 3 9 A n A \ 
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(5.10d) 


in which the nondimensionalization and the definitions of A 2 , Az and J 
have been described earlier in the inviscid full-potential solver. However, it 
should be noted that the stagnation density used for nondimensionalization is 
different from the outer-flow stagnation density. 


An approximate factorized scheme (AFl), which is a reformulation 
of an ADI scheme, is applied to solve the cavity flow. An ADI for potential 
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equation can be easily to be converted to solve the stream-function equa- 
tion. The only difference is that the density is moved from numerator to 
denominator(see Eqs. (2.1a) and (5.8a)). Consequently, there is no difficulty- 
in obtaining a converged solution for the subsonic How in the cavity. The 
numerical procedure is 

(a - 6 T1 A-j6 v )(a - 6cAi6c)C?j = (5.11) 

where C, n , = _ & ~1 /6t. and uj is a relaxation parameter. <5„(), , 

and 6s()ij are the first-order forward-difference operators defined by 

(){,; = ( )i,j + i — ( )i,y (5.12a) 

= ().>!,; -Ob; (5-126) 

and 6 r f() i ■ and 6 c() ; . are the first-order backward-difference operators de- 
fined by 

• = ( )‘J — ()t,j-i (5.13a) 

= ()t,i - (5-1-36) 


Eq. (5.11) can be derived by applying the Crank-Nicolson scheme to the two- 
dimensional heat equation. A{ and Aj are defined by 


A i = (“?) 

P J i-b,i 


A J « (“) ’ 

P'S i.j-i. 


(5.14) 


where p is updated using Eq. (5.10b). The residual is obtained from 


L K, = 


+ A 2 'H, 

T P 


+ A ?J ^ n ' 
Jp 


(5.15) 


i+k,j ' fx ,j+i 

Values of the density computed from Eq. (5.10b) are stored at cell centers, 

that is, at (i + \ ,j + ^), using values of \l/c and \J > n computed from 
(^h+iV+l ~ _ 1 + *i+ij ~ V i,j ) 


(5.16) 
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(^,) i+ jL jJ+ i = -(^i+i , j+i - * i+ i,j + #t,y+i - *ij) 


(5.17) 


Values of the density required at {i -f \,j) or (i,j 1) are obtained using 
simple averages. The ADI scheme is implemented in a two-step format: 


(« ~ 6 v A j 6 r,)n i j = auLVij (5.18a) 

where u % 2. and 

[o-TiAiT.)^, = /", ( 5 . 186 ) 


A scalar tridiagonal matrix is inverted at each step using the Thomas algo- 
rithm. A repeating sequence of a's (that is. variable time step) is used to 
speed up the convergence. A suitable sequence of a's proposed by Ballhaus 
et al. 45 is used 


<*k = a H (a L /a H ) [k l)/(M ]) k = 1, 2, 3 , ..., M 

where M is the number of elements in the sequence, for which either 6 or 8 
is used in the present study. Both an and a / v are optimized by numerical 
experimentation. 


Grid Generation 

The grids for the cavity flow is generated by an algebraic method. A 
linear function is chosen as follows, 

= S t (fl V ~' + * 6 (fl VmaX ~ 1 (5.19 a) 

y'tt,v) = yt{£)— — 37 + ybU) -- - x _ V (5.196) 

Vmax 1 V max 1 

where subscript t or 6 indicates the top or bottom wall of the cavity. The 
values of x t and yt are obtained by locally refining the outer-flow coordinates 
on the porous surface, and the values of Xf, and yi are determined by the shape 
of the cavity. 



CHAPTER VI 


NUMERICAL RESULTS AND DISCUSSION 
Inviscid Flow Solutions 

The inviscid flow solutions shown in this chapter are obtained from 
the procedure descibed in Chapter II without the entropy correction. A 14% 
thick NASA supercritical airfoil 6 and a NACA 0012 airfoil were used in a series 
of numerical studies of the effect of a porous surface on transonic performance 
with the first model. All of the computations by the inviscid-flow approach 
were performed on a 223 x 31 C-type mesh, with 162 of the 6913 grid points 
on the airfoil surface. 

Three types of porosity distribution have been examined which were 
obtained by var-ying the porosity distribution function a in Eq. (5.1). They 
are described as follows: 

Type 1 

o = constant 

Type 2 

r x ~ x \ 

V 

where x\ and 12 are the limits of the porous region shown in Fig. 1.1, and 
o max is the maximum porosity at the midpoint of the region. The Type 2 
porosity distribution is the distribution used by Savu and Trifu. 8 

Type 3 

X — X s 7T 

Xk X s 2 
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where x 3 is the horizontal position of the shock if the porous surface were 
solid, and x* represents either xi or X 2 depending on whether x is less or 
greater than x s . This function automatically adjusts the porosity distribution 
so that it decreases from the maximum value under the shock wave to zero at 
either end of the porous region. 

The performance of a 14% thick NASA supercritical airfoil with a 
porous surface has been experimentally investigated. 5-7 However, as men- 
tioned earlier, this airfoil was mounted so that only the upper half of the 
airfoil was exposed to the transonic flow. This experimental arrangement was 
modeled by reflecting the upper surface to the lower surface and computing 
the flow about the resulting symmetrical airfoil at zero angle of attack. There 
is some uncertainty concerning the airfoil coordinates in the region near the 
trailing edge where the y-coordinates of the airfoil’s upper surface are nega- 
tive. In the numerical simulation, these coordinates are modified by using a 
linear interpolation from (x,y) = (0.95, 0.0057) to the trailing edge at (1.0, 
0 . 0 ). 

To agree with the experimental shock-location, 7 the computations 
had to be performed at free-stream Mach numbers slightly lower than those 
measured. This is probably caused by the differences between the experimen- 
tal and numerical geometries including the lack of wind tunnel wall modeling 
in the computed results and the inviscid assumption. For example, the flow 
over the solid supercritical airfoil at = 0.81 in the wind tunnel can be 
approximated by using Moo = 0.795 in the computation (Fig. 6.1). Also the 
flow over the same airfoil at M 0 0 = 0.85 can be approximated by using 
= 0.805 (Fig. 6.2). Type 1 porosity with o max = 0.6 was used to model the 
experimental uniform 2.8% porosity (based on the hole area divided by the 
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Fig. 6.1 Comparison of computed and experimental pressure distribution 
on the surface of the modified NASA supercritical airfoil. 
Experimental data at M x = 0.81 and 2.8% porosity; computations 
at M co = 0.795 with d = 0.6 Type 1 porosity. 




x/c 


Fig. 6.2 Comparison of computed and experimental pressure distribution 
on the surface of the modified NASA supercritical airfoil. 
Experimental data at Moo = 0.85 and 2.8% porosity; 
computations at Moo = 0.805 with o = 0.6 Type 1 porosity. 
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total airfoil area). 3 The porous region extends from 56% to 83% in both the 
experiment and computations. The computed and experimental pressure dis- 
tributions on the porous airfoil are compared in Figs. 6. 1-6. 2. Overall, the 
agreement is quite good. The discrepancy in the trailing edge region is due, 
at least in part, to the difference in geometries and viscous effects. These 
results suggest that o maz — 0.6 may be used to simulate this flow over a range 
of Mach numbers and also indicate that the assumption of constant-cavity 
pressure results in reasonable pressure distribution predictions. 

* The variation in the pressure distribution with d is shown in Fig. 6.3. 
As a is increased, the variation in the pressure distribution becomes smaller. 
For example, as a is increased from 0.6 to 0.8, there is very little change in the 
pressure distribution. In other words, increasing the porosity beyond a « 0.6 
does hot significantly change the flow or improve the performance. 

It is difficult to make a direct comparison of the computed wave 
drag with the experimental drag presented in Ref. [7]. The measured drag data 
contain the effects of viscosity, airfoil surface roughness, and other factors that 
have not been considered in our analysis. For approximate comparison, the 
viscous drag was estimated by computing the viscous flow over the solid airfoil 
at a nominal 0.6 using the viscous/inviscid interaction code of Refs. [22- 

23]. This C | viscous ~ 0.012 was then added to the wave drag computed with 
the present code, and the results are compared with the experimental data 
in Fig. 6.4. The drag reduction effect caused by the porous surface observed 
in the laboratory is qualitatively obtained in the present study. However, the 
drag reduction occurs at lower than in the experiment. In addition, the 
increased drag at lower M ^ in the experiment is not predicted. 
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Fig. 6.3 Effects of varying the porosity strength a on the pressure 

distribution of a modified NASA supercritical airfoil at a = 0° 
and Moo — 0.795. 




L 

i 

i 

4 

L 

L 

L 

L 

L 

i 

L 

I 

c 



Fig. 6.4 Computed and experimental drag on a modified NASA 
supercritical airfoil at a = 0°. 
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The effect of a porous surface on the M = 0.8 flow past the NACA 
0012 airfoil at zero angle of attack is presented in Figs. 6.5 and 6.6. A Type 2 
porosity is used on the upper surface between x — 0.4 and x = 0.8 with o max = 
0.6. The constant-Mach number contours show a normal shock standing below 
the solid lower surface of the porous airfoil whose location and strength are not 
much different from those of the shock on the original solid airfoil. However, 
the shock above the upper surface is weakened in the presence of the porous 
surface, as shown by a group of less concentrated constant-Mach lines around 
the sonic line. Figure 6.5 also reveals that the porous surface causes a weak, 
oblique compression wave at its upstream end ( x = 0.4), due to the blowing 
from the cavity, and a readjusted compression downstream of the shock due 
to the suction of air into the cavity. The contour lines in the shock region are 
no longer normal to the airfoil surface, resulting in a lambda-shaped shock 
wave structure similar to that photographed in the laboratory (Refs. [5-7]). 

The pressure distribution on the upper airfoil surface for the case just 
presented is plotted in Fig. 6.6. The results for a porous airfoil (dashed line) 
and solid airfoil (solid line) are compared. The comparison clearly shows that 
the original steep compression through the normal shock on the solid airfoil 
has been reduced in the presence of the porous region. The original shock 
is replaced by several weaker compressions over the region covered by the 
porous surface. The resulting weaker, adverse-pressure gradient would lessen 
the possibility of flow separation. The porous upper surface has a negligible 
influence on the pressure distribution along the lower surface; that distribution 
is, therefore, omitted here. The asymmetrical pressure distribution on the 
upper and lower surfaces causes a small lift on the airfoil at zero angle of 
attack, with Cl = 0.0183. On the other hand, a 27.5% decrease in wave drag 



Fig. 6.5 Constant Mach number contours around a porous NACA 0012 
airfoil, Moo = 0.8, a = 0°, Type 2 porosity, a moI = 0.6, 
x\ = 0.4, X 2 = 0.8 (upper surface is porous). 
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from Cp = 0.0069 to 0.0050 is found for the porous airfoil. 

Much larger decreases in drag are obtained when both upper and 
lower surfaces are made porous. The effect of varying the porosity strength 
on the drag of a double-porous NACA 0012 airfoil at a = 0° is shown in 
Fig. 6.7. A Type 2 porosity distribution between x = 0.1 and x = 1.0 is used 
for all curves presented. The porous surface has a drag-reduction effect only 
when a shock appears above it at Mach numbers higher than 0.77, and that 
effect is enhanced by increasing o max (i.e., by using porous surfaces having 
smaller resistance to the penetrating flow). However, this trend diminishes 
as &ma.x reaches higher values. This is consistent with the result shown in 
Fig. 6.3. 

Changing from Type 2 to Type 3 porosity does not cause a signifi- 
cant change in drag if the normal shock appears near the center of the porous 
surface. This change assumes the same value of o max for both porosity distri- 
butions. If the shock wave is not centered in the porous region, Type 3 porosity 
is considerably better than Type 2 in smearing the shock and reducing the 
wave drag. 

The behavior of porous airfoils at a finite angle of attack is now 
described. The Mach number contours for a solid NACA 0012 airfoil at Mach 
0.75 and an angle of attack of 1° are plotted in Fig. 6.8. A shock wave appears 
only on the upper surface. A Type 3 porosity of strength d max = 0.3 is then 
distributed on 90% of the upper surface between x = 0.1 and x = 1.0. The 
resulting flow pattern and pressure distribution, plotted in Figs. 6.9 and 6.10, 
respectively, reveal that this widely distributed porosity is very effective in 
reducing the shock strength: By making the upper surface porous, Cl is 
increased from 0.240 to 0.357 while Cp is decreased from 0.0240 to 0.0008, 




Fig. 6.7 Effect of varying porosity strength on the 
0012 airfoil at a = 0°, Type 2 porosity, x 
1 2 = 1.0 (both the upper and lower surfa 
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Fig. 6.9 Constant Mach number contours around a porous NACA 0012 
airfoil, = 0.75, a = 1°, Type 3 porosity, o ma x — 0.3, 

X\ = 0.1, x 2 = 1.0 (upper surface is porous). 







which corresponds to a nearly shock-free condition. 

To study the effect of varying porosity strength on lift and drag of 
a transonic airfoil, a Type 3 porosity is distributed between z = 0.3 and x 
— 0.9 over the upper surface of a NACA 0012 airfoil. For this series of cases 
the angle of attack is fixed at 1°. The result for drag plotted in Fig. 6.11 
is analogous to that shown in Fig. 6.7 for the double-porous airfoil at zero 
angle of attack, showing that large wave-drag reductions can be achieved by 
increasing the porosity. Plotted in Fig. 6.12 are the lift coefficient data which 
indicate that lift is increased by making the upper surface porous. Unlike the 
drag coefficient, Cl is affected by the porous surface at Mach numbers less 
than 0.72 when the shock is still upstream of the porous region. The higher 
lift is caused by asymmetric changes in the pressure on the porous-upper and 
solid : lower surfaces of the airfoil. A dramatically increased lift of a porous 
airfoil was also observed by Savu et al. using the inviscid small-disturbance 
approximation (.Ref. [8]). Relatively, results from an IBL approach by Oiling 18 
were more conservative (lift was increased about 2% ). Indeed, the lift was 
improved little for a porous airfoil at lower and a, as has been indicated 
experimentally by Krogmann et al. (Ref. [9]). 

According to the present inviscid-fiow approach (although the results 
may be slightly over-optimistic and somewhat at variance with the experimen- 
tal results), it can be concluded that the shock strength can be weakened, wave 
drag can be reduced and lift may be increased by appropriate porosity distri- 
bution. In addition, the inviscid-fiow approach can provide an upper bound 
of useful levels of o of interest for viscous flow calculations. 










Viscous Flow Solution 


It is well known that viscous effects are important on airfoil per- 
formance, especially at transonic speeds. When a shock/boundary layer in- 
teraction becomes strong, the inviscid-flow analysis alone is not sufficient to 
describe the flow about an-airfoil. There are numerous examples illustrating 
the discrepancies in shock position and pressure at the trailing edge predicted 
by inviscid-flow calculations in Ref. [46]. In the present inviscid-flow approach, 
such discrepancies were noted in the previous section particularly for a solid 
airfoil. In addition, the blowing and suction at the airfoil surface have such 
a large influence on the boundary layer that v w affects the inviscid flow not 
only directly but also indirectly via the effect of S' on the v n (see Eq. (3.7)). 
The purpose of the work here is to study these viscous effects in transonic 
flow past porous airfoils. In the IBL procedure, the inviscid grid is the same 
as that mentioned earlier in this chapter, with an additional 50 grid lines in 
the rj direction for the boundary layer algorithm. In the TLNS procedure, on 
the other hand, the computations were performed on a 251 x 65 C-type mesh, 
with 163 points on the airfoil surface. 

Effects of Active Blowing and Suction 

Since the porous airfoil induces blowing in the supersonic region and 
suction in the subsonic region, it is meaningful to investigate separately ef- 
fects of blowing and suction on the shock and bound ary- layer control. The 
first test case is a NACA 0012 airfoil, at = 0.75, with a = 2.0°, the 
transition specified at 0.03 chord, and Re = 3.76 x 10 6 . Generally, the nu- 
merical results show that blowing in the supersonic region weakens the shock 
strength and smoothes the pressure gradient. However, if the blowing is too 
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strong, the results indicate that separation may occur in the blowing region 
and increase the thickness of the boundary layer approaching the trailing edge, 
which causes not only an increase in viscous pressure drag but also a decam- 
bering effect of the boundary layer leading to a decrease in lift. For the case of 
blowing ahead of the shock such as the one shown in Fig. 6.13, the results indi- 
cate that the pressure gradient is smoothed, Cp (pump drag being excluded) 
is reduced from 0.01913 to 0.01707, and Cl is also reduced from 0.37528 to 
0.30588. Even though the blowing velocity in this case is less than 1.5% of the 
free-stream velocity, C f already becomes negative in the blowing region. In 
other words, the drag reduction with normal blowing seems to be mainly due 
to the weakened shock rather than due to the boundary-layer control achieved 
by increasing the mixing rate. On the other hand, suction behind the shock 
generally increases the shock strength, moves the shock downstream and de- 
lays separation. Shown in Fig. 6.14 is an example of these phenomena. In this 
case Cl is increased from 0.37528 to 0.44354, but Cl/Cd is improved only 
slightly because wave drag and skin-friction drag are also increased. Further- 
more, the viscous pressure drag is not a dominant part of the total drag in this 
case, so the form drag cannot be reduced significantly. The boundary-layer 
control aspect is shown by the fact that C j is increased not only in the suction 
area but also in the region downstream of the suction, indicating that suction 
probably can sustain the strong shock without separation or with controllable 
separation. Therefore, in regard to drag reduction, the above numerical results 
imply that the shock can be weakened by normal blowing and boundary-layer 
control can be achieved by suction. 












Porous Airfoils 


It has been demonstrated in Refs. [22-23] that except for slight dif- 
ferences near the shock and trailing edge, good agreement can be obtained 
between results from the IBL and TLNS procedures in mildly separated tran- 
sonic solid-airfoil computations. Similarly, for a NACA 0012 airfoil at M 0 0 = 
0.8, a = 0.0° and Re — 4.09 x 10 6 , comparable Co values of 0.0125 and 0.0123 
are obtained respectively via TLNS and IBL procedures (the measured drag 
coefficient of the solid airfoil is 0.0106 47 ). Computed pressure distributions 
are plotted in Fig. 6.15(a) in comparison with measurements. The same airfoil 
is then made porous with a Type 2 porosity distribution cr max = 0.1 (the first 
model) between Xj = 0.3 and x 2 = 0.5 (on both upper and lower surfaces). 
Again comparable drag coefficients of 0.0127 (TLNS) and 0.0125 (IBL) are 
obtained numerically. Although Fig. 6.15(b) shows that the pressure jump is 
smeared by the porous surface, the pressure jump on the airfoil surface still 
has a tendency to move downstream. Since the airfoil surface over this region 
is backward facing, the lower pressure moving downstream would increase the 
pressure drag. The reduction in shock strength is not large enough to compen- 
sate the increased viscous pressure drag, so that the total drag of the porous 
airfoil becomes slightly higher than that of the solid airfoil. Such an effect 
has also been observed in the laboratory 5-7 at the lower M 0 0 range, but it 
has not been predicted by any inviscid-flow approach which can only evaluate 
wave drag. The velocity profiles in the boundary layer over the porous surface 
(from the IBL procedure) are shown in Fig. 6.16. The distance normal to the 
airfoil has been magnified ten times for clarity. It can be observed that the 
velocity profiles in the blowing region have a tendency to separate. For this 
case, the blowing (or suction) velocity at the porous surface is less than 1.0% 




Fig. 6.15 Comparison of pressure distributions for a NACA airfoil 
M oo - 0.8,o 1 - 0°, Re -• '1.00 x 10 r \ rr, - 0..'?, x- t = 0.5. 
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of freestream velocity. 

For a lifting case, a calculation (based on the TLNS procedure) was 
made with a NACA 0012 airfoil at Moo = 0.77, with a — 1.0°, Re — 4.09 x 10 6 , 
the transition fixed at 0.01, and a Type 3 porosity distribution d max = 0.07 
(the first model) from 0.378 to 1.000 which occupies a large portion on the 
upper surface of the airfoil. It can be seen in Fig. 6.17 that the pressure jump 
across the shock is weakened, Cl is increased from 0.169 to 0.183, but again 
C d is increased from 0.0125 to 0.0137 while the skin-friction drag is little 
changed. By comparison with the result for the solid airfoil, the skin friction 
is reduced in the blowing region and is increased in the suction region due to 
the porosity. 

The results in Figs. 6.18(a)-6.19(m) (obtained from the TLNS proce- 
dure) are for a symmetrical airfoil which is generated by reflecting the upper 
surface of an RAE2822 airfoil to the lower surface. The porous surfaces are 
from 0.615 to 0.805 (the second model on both upper and lower surfaces) with 
a Type 1 porosity a = 0.4. The bottom of the cavity is at y/c = 0.0. The 
flow and airfoil parameters are Re — 6.5 x 10 6 , a = 0.0° and the transition is 
fixed at 0.03. For the first case with = 0.82, Co is increased from 0.0258 
to 0.0273 in the presence of the porous surface. The effects of the porosity are 
described as follows: 

1. The comparison of C p plots(Fig. 6.18(a)) for the solid and 
porous airfoils shows that, relative to the solid airfoil, the shock is weakened 
near the airfoil surface and the pressure on the porous airfoil is lower at the 
trailing edge. The C p plots also indicate that the original one strong pressure 
jump becomes two consecutive weaker jumps, which represent the leading and 
rear legs of a lambda shock. In this case, the rear leg is slightly behind the 
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Fig. 6.17 Pressure, skin friction and normal velocity distribution on a lifting 
NACA 0012 airfoil. 
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original shock position on the solid airfoil. The lowered pressure at the trailing 
edge has also been observed experimentally by Raghunathan and Mabey who 
recently investigated the flow about a 6% thick circular arc model mounted 
on the upper wall of a wind-tunnel 11 . (See Fig. 6.18(b).) 

2. The comparison of the Mach-contour plots (Figs. 6. 18(c), (d)) 
indicates that a lambda shock structure does occur on the porous airfoil. The 
leading leg of the lambda shock slants at the beginning of the porous surface. 
This shock structure has been illustrated by the present inviscid-flow approach 
and also observed experimentally by Nagamatsu et al.(Fig. 6.18(e)) 

3. A comparison of streamline patterns near the porous region 
•’ is shown in Figs. 6. 18(f), (g). The passive flow through the porous surface 

generates a bump-like bubble (not a separation bubble) where the lambda 
shock is situated. 

4. Blowing at the leading part of the porous region reduces Cf 
and suction at the rear part of porous region makes the local C / larger in the 
suction region, as shown in Fig. 6.18(h). But downstream of the suction region, 
the flow has difficulty overcoming the adverse-pressure gradient resulting in 
values of C/ smaller than for the solid airfoil. Aft of the porous surface, Cf 
even becomes negative. 

5. Pressure contours from the cavity-flow solution are presented 
in Fig. 6.18(i). The pressure variation is very small within the cavity. Further 
description of flow in the cavity will be described in the next case. 

For the same airfoil just described but at M 0 0 = 0.85, the computed 
results indicate about a 10% reduction in drag (Cp = 0.055 for the solid airfoil, 
Cp = 0.050 for the porous airfoil). The effects of porosity are described as 
follows: 
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Fig. 6.18 Effects of porosity on a modified RAE 2822 airfoil, X\ — 0.615, 
i 2 = 0.805, Type 1 porosity on both upper and lower surfaces, 
a = 0.4, the bottom of cavity is at y/c = 0.0, Re = 6.5 x 10 , 
a = 0°, Moo = 0-82 and transition is fixed at 0.03. 

(a) Comparison of pressure distribution on the airfoil surface. 
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1. The comparison of C p plots (Figs. 6. 19(a), (b)) again shows that, relative 
to the solid airfoil, the shock is weakened near the airfoil surface and the 
pressure on the porous airfoil is lower at the trailing edge. The C p plots also 
indicate that the original one strong pressure jump becomes two weaker jumps. 
However, the second pressure jump on the airfoil surface is moved upstream 
slightly, which is opposite to the last case (M^ = 0.82 and C jr> was increased). 

The experimental results in Fig. 6.19(b) show the comparison of 
shock position between the solid and porous airfoils by Nagamatsu et al. The 
comparison of Cd versus corresponding to this figure has been shown 
in Fig. 6.4. The experimental results may indicate that the drag is reduced 
when the downstream movement of the rear leg of shock is retarded. From 
the comparison of last case (rear leg of shock moved downstream and drag 
was increased) and this case (rear leg of shock moved upstream and drag was 
reduced), there is a correlation between the experiment and the computation. 

2. Once again, the Mach-contour plots (Figs. 6. 19(c), (d)) indicate 
that the Mach contours are less concentrated around the shock on a porous 
airfoil, indicating a weakened shock. In addition, the results show that the 
position of the shock is moved upstream slightly. 

3. According to the present modeling, the effects of porosity are to 
increase the viscous losses near the airfoil and decrease the entropy changes 
away from the airfoil (referring to Figs. 6. 19(e), (f)). 

4. The separated region is enlarged on the porous airfoil (Figs. 6.19 
(g),(h)) and boundary layer becomes thicker. Fig. 6.19(h) indicates that some 
of the fluid particles that are blown out are immediately sucked back into the 
cavity. But some of the fluid particles that are blown out at the upstream 
end of the porous surface enclose a dead air region before going back into the 
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cavity. 

According to the present turbulence and porosity modeling, for the 
passive shock/ boundary layer interaction at low angles of attack, there are 
two possible flow patterns above the porous surface as shown in Fig. 6.19(i). 
There is a stagnation point in the flow field of bottom pattern (pattern 2), 
which does not occur in the top pattern (pattern l). The status of the earlier 
Fig. 6.18(g) is near a transition from pattern 1 to pattern 2, and the status 
of Fig. 6.19(h) is a typical pattern 2. Generally, when the shock strength is 
not strong , the top pattern occurs and the drag may not be reduced. When 
the shock becomes stronger, the bottom pattern may occur and the drag 
can be reduced. It can be expected that the separation bubble would burst 
intermittently and vortex shedding would occur when either Mach number or 
angle of attack is increased further, and oscillations between pattern 1 and 
pattern 2 would occur. Up to now, the only experimental investigation on 
the boundary layer near the porous region was by Krogmann, Stanewsky and 
Thiede (Refs. [9-10]). The experimental results indicate that the boundary 
layer is thickened by using a perforated surface with a cavity. 

5. As shown in Fig. 6 . 1 9 ( j ) , the suction at the rear part of the 
porous region makes C/ more negative over the suction region, which is op- 
posite to the effect shown in Fig. 6.18(h). This result should not be a surprise 
to us, in that the suction is trying to swallow the separation bubble that lies 
behind the suction region. 

6. The cavity-flow solution is presented in Figs. 6.19(k)-(m). The 
pressure variation is small within the cavity. The pressure under the shock is 
lower than at the two ends of the cavity. Since the suction area is smaller 
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x 2 = 0.805, Type 1 porosity on both upper and lower surfaces, 
o = 0.4, the bottom of cavity is at y/c = 0.0, Re = 6.5 x 10 6 , 
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Fig. 6.19(i) Streamline pattern of passive shock/boundary layer 
interaction. 
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Fig. 6.19(m) Constant Mach number contours in the cavity. 





than the blowing area, the average suction velocity is higher than the average 
blowing velocity. The highest Mach number occurs underneath the shock 
position. 

In summary, the results of the present case reveal some general effects 
at higher and a ss 0°. The shock becomes weaker, shock position is 
moved slightly upstream, and the viscous loss is increased. These findings are 
consistent with the results of the last few cases, except here the reduction in 
shock strength is more than enough to compensate for the increased viscous 
loss. Also the pattern 2 separated flow becomes mature, and downstream 
movement of the shock has been prevented. The numerical results strongly 
suggest that the total drag reduction is mainly due to a weakened shock. 
Boundary-layer separation occurs due to the blowing at the front part of the 
porous surface, and the suction at the rear part of the porous surface captures 
a dead air region. Also the present porosity model including a cavity model 
shows that the-cavity flow is creeping and that pressure variation is very small 
in the cavity. Therefore whether the first model or the second model is used 
does not affect the outer flow solution. This might not be so if very small 
cavities were considered in the interest of conserving space. 

The variation of with increasing free-stream Mach number for 
this airfoil with a slightly longer porous surface from 0.615 to 0.88 is plotted 
in Fig. 6.20. Drag reduction by the porous surface occurs at Mach numbers 
higher than 0.84, and it can be as high as 20%. Drag reductions of the same 
order of magnitude have been found experimentally by Nagamatsu, Ficarra 
and Dyer 7 in the study of a supercritical airfoil mounted on the bottom wall 
of a wind tunnel. Finally, the numerical results show that the total drag of 
the porous airfoil is reduced at higher Mach numbers while increased at lower 
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Fig. 6.20 Porosity effects on the drag of a modified RAE 2822 airfoil, 

x ] = 0.615, 1 2 = 0.88, Type 1 porosity on both upper and lower 
surfaces, a = 0.4, the bottom of cavity at y/c = 0.0, 

Re = 6.5 xlO 6 , q = 0° and transition is fixed at 0.03. 
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CHAPTER VII 


CONCLUSION 


Summary 


A numerical study is made to examine the effects of a porous surface 
on transonic airfoil performance and to help understand the flow structure 
of passive shock /boundary layer interactions. The transonic full-potential. 
1BL and TLNS algorithms modified to handle the transpiration boundary 
conditions have been used in this study. In many cases, good agreement has 
been found between experimental results and numerical results. 

In the inviscid-fiow approach, the results indicate that, by making the 
airfoil porous near the shock, the shock can be weakened to become a lambda 
shock so that the wave drag is reduced, and the lift can be increased with 
an appropriate porosity distribution. However, due partially to the neglect of 
viscous effects, the Cq versus curve, the trailing-edge pressure, and the 
shock position (without Mach-number correction) do not compare well with 
the experimental results. 

In the viscous-flow approach, computational results qualitatively ver- 
ify most of the available experimental data and improve the inviscid-flow so- 
lutions. For example, the computational results can be used to predict the 
general trend shown in the Cp versus curve, the generally lower trailing- 
edge pressure on the porous airfoils at higher Moo and a = 0°, and the thick- 
ening of the boundary layer downstream of the porous surface observed in 
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the laboratory. The results also demonstrate that the porous surface affects 
not only the shock strength and shock shape but also the position of the 
shock. Furthermore, the computational results can reveal the structure of 
passive shock /boundary-layer interactions and the streamline patterns near 
the porous surface that experiments have had difficulty with. 

Most of the work in this investigation has gone into algorithm devel- 
opment (and modification), boundary condition implementation, turbulence 
modeling, and porosity modeling (including a cavity modeling). Effort was 
also required to achieve a physical understanding of the flow patterns ob- 
tained. 

Recommendation 

In the last chapter, it has been pointed out that a porous airfoil can 
improve as well as deteriorate airfoil performance. An adaptive porosity ca- 
pability would be helpful in the effort to improve airfoil performance under 
various flight conditions! An optimal porosity distribution might be obtained 
by using a flow solver program coupled with a numerical optimization algo- 
rithm. 

Before such expectations can be realized, work needs to be done in 
turbulence modeling for massively separated flows, studying unsteady flows, 
refining the porosity model, developing a faster TLNS solver with optimum 
numerical dissipation terms, or developing a good model for the passive shock 
wave/boundary-layer interaction via the IBL procedure. Accurate experimen- 
tal data concerning the boundary layer near the porous region would be very 
helpful in constructing more realistic porosity-cavity models. 
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